load("data_for_analysis.rdata")
library(tidyverse)
library(survival)
library(rms)
library(stargazer)

#### Un-comment and run the following lines to fix stargazer error when making regression tables ####
# Source: https://gist.github.com/alexeyknorre/b0780836f4cec04d41a863a683f91b53

# detach("package:stargazer",unload=T)
# # Delete it
# remove.packages("stargazer")
# # Download the source
# download.file("https://cran.r-project.org/src/contrib/stargazer_5.2.3.tar.gz", destfile = "stargazer_5.2.3.tar.gz")
# # Unpack
# untar("stargazer_5.2.3.tar.gz")
# # Read the sourcefile with .inside.bracket fun
# stargazer_src <- readLines("stargazer/R/stargazer-internal.R")
# # Move the length check 5 lines up so it precedes is.na(.)
# stargazer_src[1990] <- stargazer_src[1995]
# stargazer_src[1995] <- ""
# # Save back
# writeLines(stargazer_src, con="stargazer/R/stargazer-internal.R")
# # Compile and install the patched package
# install.packages("stargazer", repos = NULL, type="source")
# library(stargazer)

#### Figure 2 ####
ggplot(d_state_year_FSA_ARCHV,
       aes(x = year, y = parse_number(newsletter_phase))) +
  geom_line(size = 1) +
  scale_x_continuous(breaks = seq(min(d_state_year_FSA_ARCHV$year), max(d_state_year_FSA_ARCHV$year), 3), name = element_blank()) +
  scale_y_continuous(breaks = 0:3, labels = c("Phase 0", "Phase 1", "Phase 2", "Phase 3"),
                     name = element_blank()) +
  geom_segment(aes(x = year, xend = year,
                   y = -Inf, yend = Inf,
                   color = SED_type == "1_elective_aspirant"),
               size = 3, alpha = .3) +
  scale_color_manual(name = "Type of incumbent SED", labels = c("Other type", "Elective aspirant"),
                     values = c("gray80", "black")) +
  facet_wrap(~ str_to_title(state), ncol = 5) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust = 0.5),
        legend.position = "bottom")

#### Figure 3 ####
ggplot(d_state_year_FSA_ARCHV %>%
         group_by(year) %>%
         summarise(n_states_phase_3 = sum(newsletter_phase == "phase_3", na.rm = T),
                   n_states_phase_2_or_3 = sum(newsletter_phase %in% c("phase_2", "phase_3"), na.rm = T)) %>%
         ungroup()) +
  geom_line(aes(x = year, y = n_states_phase_3),
            linetype = "solid", linewidth = 1) +
  geom_line(aes(x = year, y = n_states_phase_2_or_3),
            linetype = "dashed", linewidth = 1) +
  scale_x_continuous(breaks = min(d_state_year_FSA_ARCHV$year):max(d_state_year_FSA_ARCHV$year), name = element_blank()) +
  scale_y_continuous(limits = c(0, 50),
                     breaks = seq(0, 50, 2),
                     name = "Number of states issuing state FSA newsletters") +
  annotate(geom = "text", x = 2005.5, y = 3.5, label = "Nationally standardized state newsletters\n(Phase 3)", angle = 7.5, size = 3) +
  annotate(geom = "text", x = 2003, y = 10.5, label = "Any state newsletters\n(Phase 2 or 3)", angle = 10, size = 3) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust = 0.5),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(size = .05, color = "black"))

#### Table 1 ####
msurv_1.1 <- coxph(Surv(month_index, reach_phase_3) ~ SED_type, 
                   cluster = state,
                   d_surv_month_ARCHV)
msurv_1.2 <- coxph(Surv(month_index, reach_phase_3) ~ SED_type + age + previous_SED_appointment, 
                   cluster = state,
                   d_surv_month_ARCHV)
msurv_1.3 <- coxph(Surv(month_index, reach_phase_3) ~ SED_type + age + previous_SED_appointment + 
                     I(n_bordering_states_phase_3 / n_bordering_states * 100), 
                   cluster = state,
                   d_surv_month_ARCHV)
msurv_1.4 <- coxph(Surv(month_index, reach_phase_3) ~ SED_type + age + previous_SED_appointment + 
                     I(n_bordering_states_phase_3 / n_bordering_states * 100) +
                     I(census_principal_operators_whose_primary_occupation_is_farming / pop * 100) +
                     I(census_farms / pop * 1000),
                   cluster = state,
                   d_surv_month_ARCHV)
msurv_1.5 <- coxph(Surv(month_index, reach_phase_3) ~ SED_type + age + previous_SED_appointment + 
                     I(n_bordering_states_phase_3 / n_bordering_states * 100) +
                     I(census_principal_operators_whose_primary_occupation_is_farming / pop * 100) +
                     I(census_farms / pop * 1000) +
                     internet_access + n_secretary_disasters, 
                   cluster = state,
                   d_surv_month_ARCHV)

stargazer(msurv_1.1, msurv_1.2, msurv_1.3, msurv_1.4, msurv_1.5,
          dep.var.labels = "Adoption of standardized state FSA newsletters",
          covariate.labels = c("Elective aspirant", "Ambitious bureaucrat", "Revolving door lobbyist", "Experienced politician", "Policy expert",
                               "Age", "Previous SED appointment",
                               "Adoption rate among bordering states",
                               "Percentage of population in farming", "Number of farms per thousand population", 
                               "Internet access", "Number of declared natural disasters"),
          apply.coef = exp, p.auto = FALSE, t.auto = FALSE, report=('vc*p'), digits = 2,
          df = F, star.cutoffs = .05, notes.append = F, notes = "Standard errors clustered by state. *$p<0.05$",
          type = "text")

#### Table 2 ####
m_full_text_1 <- ols(participation_encouraging_words ~ SED_type +
                       count_words +
                       state + year,
                     d_state_year_letters_ARCHV,
                     x = T, y = T) %>%
  robcov(cluster = d_state_year_letters_ARCHV %>%
           pull(state))
m_full_text_2 <- ols(participation_encouraging_words ~ SED_type +
                       count_words + 
                       I(census_principal_operators_whose_primary_occupation_is_farming / pop * 100) +
                       I(census_farms / pop * 1000) +
                       n_months_secretary_disaster +
                       state + year,
                     d_state_year_letters_ARCHV,
                     x = T, y = T) %>%
  robcov(cluster = d_state_year_letters_ARCHV %>%
           pull(state))
m_full_text_3 <- ols(participation_encouraging_words ~ SED_type +
                       count_words + 
                       I(census_principal_operators_whose_primary_occupation_is_farming / pop * 100) +
                       I(census_farms / pop * 1000) +
                       n_months_secretary_disaster +
                       state + party_of_pres + party_regime,
                     d_state_year_letters_ARCHV,
                     x = T, y = T) %>%
  robcov(cluster = d_state_year_letters_ARCHV %>%
           pull(state))

m_all_programs_1 <- ols(participation_encouraging_program_references ~ SED_type +
                          count_program_references + 
                          state + year,
                        d_state_year_letters_ARCHV,
                        x = T, y = T) %>%
  robcov(cluster = d_state_year_letters_ARCHV %>%
           pull(state))
m_all_programs_2 <- ols(participation_encouraging_program_references ~ SED_type +
                          count_program_references + 
                          I(census_principal_operators_whose_primary_occupation_is_farming / pop * 100) +
                          I(census_farms / pop * 1000) +
                          n_months_secretary_disaster +
                          state + year,
                        d_state_year_letters_ARCHV,
                        x = T, y = T) %>%
  robcov(cluster = d_state_year_letters_ARCHV %>%
           pull(state))
m_all_programs_3 <- ols(participation_encouraging_program_references ~ SED_type +
                          count_program_references + 
                          I(census_principal_operators_whose_primary_occupation_is_farming / pop * 100) +
                          I(census_farms / pop * 1000) +
                          n_months_secretary_disaster +
                          state + party_of_pres + party_regime,
                        d_state_year_letters_ARCHV,
                        x = T, y = T) %>%
  robcov(cluster = d_state_year_letters_ARCHV %>%
           pull(state))

stargazer(m_full_text_1, m_full_text_2, m_full_text_3, m_all_programs_1, m_all_programs_2, m_all_programs_3,
          dep.var.labels = c("Number of participation-encouraging words", "Number of participation-encouraging references"),
          covariate.labels = c("Elective aspirant", "Ambitious bureaucrat", "Revolving door lobbyist", "Experienced politician", "Policy expert",
                               "Total number of words in mega-newsletter",
                               "Percentage of population in farming", "Number of farms per thousand population", 
                               "Number of months with declared natural disaster",
                               "Total number of program references in mega-newsletter",
                               "Republican president", "Divided party control of Congress", "Opposition party control of Congress"),
          add.lines = list(c("State FEs", rep("Y", 6)),
                           c("Year FEs", c("Y", "Y", "N") %>% rep(2))),
          omit = "state.*|year.*", df = F, digits = 2,
          star.cutoffs = .05, notes.append = F, notes = "Standard errors clustered by state. *$p<0.05$",
          type = "text")